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ABSTRACT 


The complex geometry of many parts of a ceramic gas tur- 
bine engine and the need to know the stresses in these parts 
to a high degree of accuracy indicate that the finite element 
method should be employed. One analysis code used a finite 
element program based on an isoparametric 12-node brick. 

This allowed quadratic variation in one coordinate direction, 
but only linear variation in the other two coordinate direc- 
tions. This introduced an artificial stiffness into the 
structure. This study was made to determine the effects of 
using a totally quadratic isoparametric brick, vice one 
which is quadratic in only one direction. 

A distributed load preprocessor was also developed and 


tested. 


FORTRAN IV was used throughout. 
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I. INTRODUCTION 


A. WHY CERAMICS? 

There are three major advantages in using ceramics in 
the hot section of a gas turbine: improved engine perform- 
ance, higher durability, and less dependence on expensive 
imported materials. 

Simple Brayton cycle analysis shows that a higher tur- 
bine inlet temperature (TIT) results in a higher thermal 
efficiency, and lower specific fuel consumption. Blade 
material considerations usually limit TIT. Intricate blade 
cooling schemes, such as those used in the LM2500 and the 
FT9, use impingement and convective heat transfer and exter- 
nal film cooling, to lower the temperature of the turbine 
blade. This degrades engine performance and introduces 
manufacturing complexities. Additionally, the cooling 
holes are susceptible to plugging, which may lead to even- 
tual failure. 

. Ceramics, however, can withstand higher turbine inlet 
temperatures without blade cooling. Ceramic components 
also have significantly higher resistance to hot corrosion 
than metal "superalloy" engines. The corrosion properties 
of ceramics make them better candidates to withstand the 
higher fuel impurities expected in the future. 


At present, metal engines require expensive, imported 


metals such as chromium, cobalt, columbium, and nickel. 


These metals are in limited supply and come from areas of 
potential political instability. Ceramics rely on inexpen- 


sive raw materials which are found in the United States, 


B. DESIGN CONSIDERATIONS 

Ceramics are brittle materials. This does not mean 
that they are necessarily fragile. A brittle material will 
fracture with little or no plastic deformation. <A fragile 
material will tracture at low stress with little or no 
plastic deformation. It has been suggested that the word 
"inductile"” be used instead of brittle (1). The inductile 
nature of ceramics is the crux of the design problem with 
this material. Local stress concentrators, such as the tip 
of a crack, will not vield in plastic deformation. If the 
stress concentration is not relieved, material fracture may 
result. 

Another major problem is the variability of data on the 
material properties of ceramics, 

These two problems have a major impact on the design. 
Design procedures that were valid for most engineering 


metals may fail when applied to ceramics. 


Cy. PROJECT BACKGROUND 

Some of the calculations performed later in this work 
are concerned with the design of a ceramic engine tor a 
Defense Advanced Research Projects Agency CDARPA) project. 


The prime contractor, the AiResearch Division of the 


Garrett Corporation, used their T76 model metal engine as 
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the basis for the ceramic design. They indicated that by 

replacing certain metal hot section parts, they could im- 

prove engine output power by 40%, and reduce specific fuel 
consumption by 10%, as shown in Fig. 1. An engine cross- 

section (Fig. 2) shows the ceramic hot section. 

As mentioned earlier, the brittle nature of ceramics 
drives the design. In the case of the Garrett engine, a 
compliant layer was added in the contact region of the 
dovetail as an attempt to reduce the contact stresses. 
Blade tip clearance was increased to reduce the possibility 
of impact of the blade on the shroud. In addition, a prob- 
abilistic approach, based on the Weibull distribution, 
was used as part of a design methodology that was aimed at 
predicting probability of success. These are but a few of 


the examples of how ceramics affected the design. 


D. REASON FOR THESIS 

The complex geometry of many of the ceramic engine's 
parts, and the need to know the stresses to a high degree 
of accuracy, indicated that the finite element method 
should be employed. One analysis code used a finite element 
program based on an isoparametric 12-node brick. This 
allowed quadratic variation in one coordinate direction, 
but only linear variation in the other two coordinate direc- 
tions. This introduced an artificial stiffness into the 
structure. This study was made to determine the effects of 
using a totally quadratic isoparametric brick, vice one 


which is quadratic in only one direction, 
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IZ. ‘THEORY 
Recall from strength of materials: 
M = EI/p 


where M = moment, E = Young's modulus, I = moment of iner- 
tia, and p = radius of curvature. 

Consider a beam made of an isotropic material, with a 
constant cross~section. If the beam is in pure bending, 
the moment is constant across the span. If Young's modulus 
and moment of inertia are also constant, the radius of cur- 
vature is a constant. A curve which has a, constant radius 
of curvature is an arc of a circle. 

Figure 3(a) shows a 4-node plane element. There can 
only be linear variation on each of the sides. If the ele- 
ment is deformed in pure bending, each of the sides remains 
straight. The top and bottom surfaces do not reflect any 
of the curvature expected. This is a very stiff approxima- 
tion. 

A 6-node plane element, shown in Figs. 3(b) and 3(c), 
yields a better approximation to the bending phenomenon, 
depending on the orientation of the one quadratic side. 
Orientation 1 does not reflect any of the bending curvature. 
Orientation 2, however, does. The finite element method 
would represent the deflected curve, which is defined by 3 


nodes, by passing a parabola through the node points. This 


15 


can give a very close approximation to the circular arc 
predicted by strength of materials. It should be noted, 
however, that orientation 2 is constrained by the same 
"planes remain plane" restriction applied to slender beam 
theory. Thus, depending on the geometry and loading, a 6- 
node plane element with orientation 2 may or may not accu- 
rately reflect bending for a deep beam. 

A totally quadratic element, shown in Fig. 3(d), is 
not constrained by orientation, geometry or loading. An 8- 
node plane element can give a better representation of bend- 
ing. 

For a given material, a stiffer element, such as the 4- 
or 6-node plane element, will deform less than it should. 
Since strain is based on displacement, and stress is based 
on strain, the stiffer structure will predict a lower stress 
than is actually present. This is a potentially dangerous 
design situation. 

The development presented here is based on a two-dimen- 
sional plane element. Similar reasoning can be applied to 
three-dimensional bricks. A 12-node brick, for example, 
has only one quadratic Side. Depending on the loading, 
geometry, and orientation of the quadratic side, a 12-node 
brick could approach the correct solution. It will always, 
however, be stiffer than the actual structure. If it is 
oriented incorrectly with respect to the applied moment, it 


could result in seriously incorrect answers for a coarse 


mesh system. In addition, since two of the coordinate 
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directions are stiffer than the third quadratic coordinate 
direction, an artificial anisotropy is introduced into the 
problem. 

A fully quadratic brick does not have the problems men- 
tioned above. Although the discretization process intro- 
duces some artificial stiffness, it is very slight. A 
fully cubic brick would even be a better representation of 


the actual structure. 
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III. BENDING SPECIMEN PROBLEM 


To evaluate the accuracy of a finite element computer 
program, it is often useful to compare its results to experi- 
mental results for a problem that can be verified experimen- 
tally. A 4-point flexure test has been used very often for 
this purpose [3]. This chapter explains how the flexure 
test can be modeled using 12-node and 20-node bricks. The 


results are also discussed. 


A. FOUR POINT FLEXURE TESTS 

Figure 4 shows a schematic of a 4-point flexure rig. 
It is used to provide test results from a specimen in pure 
bending. Figure (5a) shows a free body diagram for the 
specimen. Figures (5b) and (5c) show the shear and moment 
diagrams. The moment diagram shows that the portion of the 
Span between the applied forces on the top surface has no 


shear, and maximum constant moment, i.e., pure bending. 


B. MODELING TECHNIQUES 

A standard 0.125 inch x 0.250 inch x 3.0 inch bend 
specimen was used. The specimen was modeled using a vari- 
ably-noded three-dimension brick from the ADINA [4] and 
SAP IV [5] finite element programs. Maximum use was made 
of the symmetry in the problem. Only 25% of the specimen 
needed to be modeled, as shown in Figs. 6a,b,c. Symmetry 


in the spanwise and long tranverse directions was used. 


A cross-sectional view of the deformed bar, Fig. 7, shows 
the symmetry in the long transverse direction. The mid-span 
cross-section was given a guided end condition (no displace- 
ment in the x-direction). The midplane in the long trans- 
verse direction was also given a guided end condition (no 
displacement in the y-direction). 

Since the model was very regular in shape, automatic 
nodal point coordinate and connectivity generation were ex- 
tensively used. Both were generated in the spanwise direc- 
tion to take maximum advantage of the automatic generation 


features of the codes employed. 


C. FLEXURE CASE STUDY 

A 12=node brick has one quadratic direction. The other 
two directions are linear. When used to model the bending 
specimen, three distinct orientations of the brick can occur. 
Figure 8 shows each of these orientations with respect to 
the bending moment. An analysis was performed using the 
SAP IV program for each of the three orientations, called 
Case 1, Case 2, and Case 3. The applied forces were modeled 
as distributed loads. 

1. Analysis of Results 

The three-dimensional finite elements used to solve 

the problem would give a solution converging to an exact 
three-dimensional elasticity theory solution. However, such 
4 solution in closed form does not exist. Therefore, the 


results were compared to the usual strength of materials 


beam theory. Such a theory is not perfect and circumspection 
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had to be used in such a comparison. For instance, beam 
theory does not predict any anticlastic curvature [18], but 
the finite element solution shows some. However, due to 
the boundary conditions and symmetry, beam theory should 
apply to the x-z plane of symmetry with a very good accura- 
ey. The displacements of the neutral axis of this plane 
from x = 0 to x = L were plotted versus the span length. 
Displacements were normalized by dividing them by the mid- 
span (x = 0) displacement predicted by beam theory (5,)- 
This facilitates comparisons in terms of the percentage of 
the maximum deflection of beam theory. 

Two loadings and materials were used. The first 
specimen was made of steel, and had a 300 1b load. The 
second specimen was hot pressed silicon nitride (HPSN), 
and had a 180 lb load. The mechanical properties of HPSN 
are given in Ref. [6]. 


The deflection curves were calculated using [7]: 


A w<xea>® 


o 4 o, + se * eee 


r w(L-a) 2 2 
where Sq + = (2L° + 2aL - a’) 
Ma = w(L-a) = reaction end moment 


w = applied force 


= Young's modulus 


os 
' 


= moment of inertia 


Figure 9 shows that L is half the distance between supports 


on the lower side of the specimen, and w is the total load 


applied, i.e., 300 1b or 180 lb. The symbol < > denotes a 


ee 


discontinuous function, i.e., for x<a the quantity <x-a> = 0, 
for x>0 <x-a> = x-a. 

The results of these analyses are shown in Figs. 10 
and ll. The Case 1, 2 and 3 curves correspond to three 
orientations of the 12-node brick. Case 1, the orientation 
with the quadratic side along the axis of bending, is the 
closest of the three to the values predicted by beam theory. 
The maximum deflection of the ceramic specimen is approxi- 
mately 96% of the deflection predicted by beam theory. This 
model is slightly stiffer than beam theory. Cases 2 and 3, 
which have the quadratic side perpendicular to the axis of 
bending, are significantly stiffer models than Case 1. The 
maximum deflection predicted by these cases is only 70% of 
the deflection predicted by beam theory. 

In the displacement-based finite element method, 


stress and displacement are related by [8]: 


{e} (B) {u} 


and {t} = (C}{e} = [(C) (B){u} 


where it) a vector composed of the six independent 

components of the stress tensor. 

{e} = a vector composed of the six independent 
components of the small strain tensor. 

{B]) = strain-displacement transformation matrix. 


{u} = displacement matrix. 


{C)} = matrix of material properties. 


This formulation clearly shows that a stiffer struc- 
ture (i.e., less displacement and displacement gradients) 
will produce less stress. 

Case 4, shows results which are slightly more flex- 
ible than beam theory. Since beam theory is inherently 
stiff due to the approximations made during its development 
(i.e., negligible shear stress), the Case 4 model is closer 
to predicting the actual stress distribution than either 
Cases l, 2, 3 or beam theory. It should also be remembered 
that the actual structure modeled will bend in all three 
directions. Clearly, the only way to fully represent such 
a deformation is to use a fully quadratic element. Any 


element that is less than fully quadratic will be too stiff. 


D. TEMPERATURE GRADIENT CASE STUDIES 

The temperature gradient case studies used the same 
models as the flexure case study. Cases 1, 2, 3 refer to 
the same orientation of the 12-node brick as in the flexure 
study. Case 4 is the 20-node brick model. 

Hot pressed silicon nitride was used as the specimen 
material. The coefficient of thermal expansion is given in 
Ref. (6). A 100°F temperature gradient was imposed across 
the short transverse (0.125 inch). The stress-free refer- 
ence temperature was 1SOO°F. A variably-noded (12~21) 
three-dimensional brick from the ADINA finite element pro- 
gram was used to model the problem. Nodal temperatures were 
input via a temperature tape. This procedure is described 


in Appendix A, 
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As in the flexure case study, the normalized deflections 
of the y-z plane of symmetry were plotted versus distance 
along the span. Beam theory was used as the basis for com- 
parison. The deflection predicted by beam theory is [7]: 


= a 2 
6 = Sa + OE (T, - T,) x 


Re 2 
where: 6, = -aL (T, ~ T,)/2t 
a = coefficient of thermal expansion (inch/inch/°F) 
t = depth of beam 
(T, - T,) = temperature difference through depth of beam 


L = half the distance between the lower reactions 


Two different case studies evolved from the temperature 
problem. Reference [9] states that the Bernoulli-Euler as- 
sumptions can be used in practical analyses of beams under 
thermal loadings. The assumptions state, in part, that the 
effects of lateral contraction can be ignored, i.e., Poisson's 
ratio (v) may be taken as equal to zero. In the case studies 
that follow, Poisson's ratio was used in the first study, and 
set equal to zero in the second. 

1. Results of the Case Study using Poisson's ration (v) 

The normalized deflection curves for this case study 
are shown in Fig. 12. The displacements from x=0 to x=L are 
plotted. Cases 2 and 3 are stiffer than beam theory. The 
curves show that the maximum deflection calculated for these 
two cases is 75% of the maximum deflection predicted by beam 


theory. 


The results of Case 4, which used the fully quadratic 
20-node brick, show excellent agreement with simple theory. 

Figure 12 shows that the results of Case 1 (the 12- 
node brick with quadratic side along the span) predict a de- 
flection greater than theory. It is puzzling that the re- 
sults show that the element is even more flexible than the 
20-node brick. The temperature algorithm used by ADINA was 
checked extensively. No errors were detected. No satis- 
factory explanation for the 12-node brick behavior has been 
found. It is suspected that the inherent anisotropy of the 
l2-node brick may be the cause of the problem. However, 
extensive further examination would be necessary to resolve 
the dilemna. 

Possibly more interesting than the displacement re- 
Sults are the stresses calculated. Recall that if a linear 
temperature gradient were imposed on a bar that had no ex- 
ternal constraints, the bar would deform into an are of a 
circle. Also there would be no stress in the bar according 
to a fundamental theorem of elasticity. In this case study, 
the bottom nodes at x = 0.75 inch were fixed in the z-direc- 
tion. The part of the specimen to the left of this con- 
straint (i.e., toward the free end) should, however, show no 
stress until one reaches the immediate vicinity of the con- 
straint. Each of the cases that use the 12=-node brick pre- 
dict stresses, both normal and shear, on the order of 10° 
to 10" PSI in the free end. The 20-node brick, however, 
predicts stresses of approximately 5 PSI or less in this 


: 6 : - = Do whe ge > 
region. (£245x10° psi, v0.2, azl.?xl0° in/in/°F). 
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The above discussion concerning the 12-node brick 


models indicates the possible uncertainty that can result 
from using this element in thermal calculations. 
2. Results of Case Study with v = 0 

This case study was undertaken in an attempt to ex-~ 
plain the results of the Case 1 orientation of the 12-node 
brick discussed above. Everything was the same as in the 
previous study, except that the effects of Poisson's ratio 
were neglected. The results are shown in Fig. 13. The Case 
2 and 3 models have become significantly stiffer. The maxi- 
mum deflection was approximately 67% of the value predicted 
by theory, as. compared to approximately 75% when Poisson's 
ratio was used. The results for Case 4 have remained the 
same (to three significant digits). The displacements for 
Case 1 have decreased to such a degree that they were almost 
identical to Case 4. The spurious free end stresses cited 
previously remained. 

In summary, by ignoring the effects of Poisson's 
ratio, three of the models have become significantly stiffer, 
while the fully quadratic model has remained substantially 
the same and in close agreement with the theory. 

It is interesting to note that if the specimen were 
truly unconstrained and had a linear temperature gradient, 
the results of analyses that take into account Poisson's 
ratio should agree exactly with the results of analyses that 


do not [9]. 


IV. AIRFOIL MODELING 


This chapter discusses the techniques developed to con-~ 
struct a finite element mesh for the airfoil section of a 
ceramic gas turbine blade. This particular engine airfoil 
lent itself to detailed analysis because of the relative 
ease with which the mesh could be generated. This is due 
to the linear fairing between the tip and root sections of 
the airfoil. Although linear fairing was not optimal from 
an aerodynamic point of view, it was desirable from a manu- 
facturing point of view. Its relative ease of fabrication 


caused it to be used in practice. 


A. DEVELOPMENT OF AIRFOIL PROFILES 

The Garrett Corporation provided Mylar cross-sectional 
profiles of the first stage rotor blade of their experimen- 
tal ceramic turbine [10]. These profiles were for various 
radial distances as measured from the axis of rotation 
along the stacking axis of the blade. In addition to an 
outline of the cross-section and its relation to the stack- 
ing axis, a table of x-y coordinates was given for the pro- 
file. These coordinates described 62 points on the profile. 

The 62 profile points were assembled into x and y ma- 
trices. These matrices were used to develop cross-sectional 
plots on the time sharing terminal and later on the Calcomp 


plotter. Figures 14 and 15 show representative cross- 
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sections of the tip and hub sections. The profile points 
are marked by an 'x' and are connected by straight line 
segments. 

Since there was linear fairing between the tip and hub, 
it could be easily verified that the coordinates listed on 
the Mylar drawings were given at the same relative positions 
along the foil profile. In other words, the ith point on 
the R=3.1 inch profile corresponded to the ith point on the 
R=4.0 inch profile, where R is the radial distance along the 
stacking axis from the axis of rotation. 

The root profile was then divided into elements. Care 
was taken to develop elements with convex corners. Also, a 
relatively fine mesh was desirable in the leading and trail- 
ing edge areas of the airfoil to enable accurate modeling 
of the anticipated stress concentrations there. The result~ 
ing leading and trailing edge elements proved quite trouble- 
some; this will be discussed more fully below. 

The mesh developed in the airfoil root section dictated 
the mesh in the tip section. This was due to the unique 
correspondence of points in the tip and hub sections. The 
tip section was then connected to the hub section by linear 
nodal coordinate generation. Connectivities were generated 
from the leading to the trailing edge. This scheme took 
maximum advantage of the automatic generation features of 


SAP IV and ADINA. 


B. THE NEGATIVE JACOBIAN DETERMINANT 
As mentioned earlier, the leading and trailing edge ele- 


ments were ill-behaved. The leading and trailing edges of 


this blade are arcs of a circle. When an element with qua- 


dratic sides is used, these elements are prone to re-entrant 
angles, and resultant negative or zero determinants of the 


Jacobian matrix. In other words, the inverse of the Jacobi- 


an matrix does not exist because there is not a unique cor- 


respondence between the natural and local coordinates of the 


element [8]. This will stop the program because the inverse 
of the Jacobian matrix is used in the strain-displacement 


Matrix [B]. Recall: 


{e} = (B){u} 


a . 
— ) = (7 | 
%4 
oe: 
a2 


2) x: 
J gl 


\ T 
Also: {K] = Svon (B] (C] {BJ} d¢VOL) 


[J]~~ = inverse of the Jacobian matrix 
[K] = stiffness matrix 

{u} = (HI {u,} - displacement field 
{H] = matrix of shape functions 

[C] = matrix of elastic properties 
{e} = strain matrix 


r,s,t = natural coordinates 
As the above equations show, the inverse of the Jacobian 
is an essential part of the calculations for the strain ma- 
trix, and the stiffness matrix. Experience has also shown 
that when the determinant of the Jacobian is small (107? or 


smaller) the displacements at that node may lose almost all 


their numerical significance. 


C. MESH VERIFICATION THROUGH PSAP1 

Kibler [11] implemented a finite element preprocessor 
for mesh verification for the SAP IV and ADINA programs, 
This enables the user to check visually the coordinates and 
geometry of a mesh by generating graphical displays of the 
completed mesh. Oblique orthographic projections are pro- 
duced using the NPS Calcomp model 765 plotter. Other dis- 
plays are possible. 

The PSAP1 program is especially useful because of the 
various mesh presentation options in the program. Among 
them are: 

1. Node numbering 


2. Element numbering 
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3. Exploded plots 

4, Full rotation about each axis 

5. Sectioned plots based on geometry, 
node number, or element number 

6. Displacement postprocessing using either 
the deformed structure, or displacement 
vectors at the nodes. (See Chapter VII) 

Other options, as well as a User's Manual, are fully dis- 
cussed in Ref. [ll]. 

Figures 16 through 19 show the various mesh presentations, 
for the linear airfoil sections. PSAP1 was used extensively 
in the production of the airfoil mesh. The clear graphical 
presentation assures the engineer that he is using a mesh 
free of geometric or connectivity errors. 

The elements with negative Jacobians were separately 
selected from the mesh and greatly magnified. The re-entrant 
angles were clearly visible on these plots as shown in Fig. 
20. Large scale plots of the hub and tip cross-sections for 
these elenente were made by hand. The mid-size nodes were 
moved inward toward the geometrical center of these elements. 

Two airfoil meshes were developed. One was composed of 
20-node bricks, and the other of 12-node bricks. The 12-node 
bricks were oriented such that the quadratic side could be 
used to model the airfoil camber. This orientation was the 
same as the one used by Garrett. Since the 12-node brick 
has only one quadratic side, part of the leading and trailing 


edge curves were eliminated. This is shown in Fig. 21. Aside 
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from this difference, the geometry represented by the two 


meshes was identical. The 12-node airfoil is shown in Fig. 


22. 


The characteristics of each of these meshes are shown 


in Table I. 
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V. AIRFOIL CASE STUDY 


After thoroughly investigating the ceramic bar problem, 
the next step was to analyze linear airfoil section of the 
ceramic gas turbine blade using the ADINA code 12-node and 
20-node bricks. The bar problem showed that significant 
differences occurred as a result of the brick type and 
orientation used. The purpose of the airfoil study was to 
determine the difference, if any, between the 12-node brick 
airfoil and a 20-node brick airfoil. Each airfoil was 
loaded using the actual design point conditions that could 
be expected in the ceramic gas turbine. Two general types 


of loading were investigated: centrifugal and thermal. 


A. CENTRIFUGAL LOAD 

A centrifugal load preprocessor was developed at NPS, 
and is discussed fully in Ref. [12]. The preprocessor can 
be used to generate consistent nodal loads for a given 
angular speed. The rotor design speed of 41,730 RPM was 
used. 

Figures 23, 24 show the deformed foil superposed on the 
undeformed foil. (Chapter VII discusses using PSAP1 to plot 
the deformed structure.) The root plane coincides exactly 
in all cases because the nodes there were fully constrained. 
Two general deformations can be seen: first, there is a 


general elongation of the blade in the radial direction, and 


32 


oe et RRR ——— ee 


second, there is an untwisting of the blade. Figures 25-26 
show the deformed 12-node and 20-node foils superposed on 
each other. It can be seen that the 12-node model is stif- 
fer in representing the untwisting deformation. Tables II 
and III show comparisons of stresses and maximum principal 


stresses, as calculated by the 12-node and 20-node bricks, 


at selected nodes in the root section. At each of the tabu- 
lated nodes, the stresses predicted by the 20-node brick are 
higher than those predicted by the 12-node brick. The 
stress contours for the root section are shown in Figs. 29 


and 30. 


B. THERMAL LOAD 

Reference [13] shows the steady state isotherms. For 
this study, the isotherms were modeled as constant for 
various radial distances along the foil. The cited reference 
shows that this is a good assumption in the root section, 
where the isotherms are flat and closely spaced. In the 
upper portion of the foil, however, the isotherms have a 
definite curvature in the radial direction. Figure 31 shows 
the nodal temperatures used in this study. Preparation of 
the nodal temperature tape is discussed fully in Appendix A. 

Figure 32 shows the deformed foil superimposed on the 
original geometry. A general enlargement of the foil in all 
directions is seen. Table IV shows a comparison of thermal 
stresses at selected points. (It should be noted that these 


are not evaluated at nodal points, but rather at the Gauss 


points. ADINA does not give nodal stresses for thermal load- 


te i ain nee Nn haa eee ew. 
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C. CONCLUSIONS FROM AIRFOIL STUDY 
This comparative study shows how the choice of a finite 


element can affect the results. The 12-node brick, even 


though oriented in the most advantageous way, is stiffer 


than the fully quadratic three-dimensional element. The 
stiffer element produces lower stresses in the root section. 
It is these root section stresses which transmit the tensile 
forces and moments to the blade platform and eventually to 


the critical dovetail area. 


VI. DISTRIBUTED LOAD PREPROCESSOR 


The version of ADINA installed at NPS does not contain 
a distributed load feature. Finite element texts, such as 
Ref. (14), often state a consistent allocation of a uniform 
surface load for a rectangular element. The consistent 
load for a curved surface, such as the airfoil studied in 
this report, is more complex. Since a distributed load is 
used to model the aerodynamic pressure on an airfoil, a 
program to calculate consistent pressure loads for curved 


surfaces was neces sary. 


A. DEVELOPMENT OF THEORY 
A consistent load is defined as a set of nodal forces 
which produce the same amount of work as the original force. 


Stated mathematically: 
2 2 2S 


where; <v;> = vector of nodal consistent loads 
{u;} = nodal displacements 
<f*> = vector of surface forces 
{u } =(Ni} tu} = displacement field 


IN;] = shape functions 
Consider the general surface shown in Fig. 33. 


Work = J, n pdAsu = <v.> ae (1) 


A 
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unit normal to the surface (positive outward). 


pressure on dA. 


displacement vector = ui + vj + wk 


ease of computation, it is desirable to perform the 


integration in the natural coordinate system. The vector 


ndA is transformed to the natural coordinate system by: 


dR dR 
ndA = Gee x 38) da dg 


where: 


i) 
“ 


xi + yj + zk = position vector of point on 
surface as shown in Fig. 33. 
da d&g = differential area in natural coordinate 
system. 
a,8 are dummy variables standing for the appropri- 
ate natural coordinates which define the 


pressure face (See Table V). 


Q 
In 


This is a vector tangent to a line a where 8 


2 


is constant. The components of this vector 


are elements of the Jacobian matrix. 


aR 
B is similarly defined. 


The pressure and displacement fields can be discretized 


with the same shape functions: 


Pp = <N;>{p;} (2) 
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u N; 0 0 us 
Vv = 0 N; 0 Vz C3) 
W 0 0 N; Ws 
where: H, * shape functions 
{p,} = matrix of nodal pressures 
= matrix of nodal displacements 
Let 
aR aR 
(qe x Gq) da dB = CAVE + AQ] + AK) da dB (4) 
then: 
ndAsu = <A, ho A3> Vv da d6 
Ay * ®o 4 “9.9 ~ "2,0 oR 
Ao ™ Bag “1.8 “ie “i568 
Ag * %. 4 %2,8 ~ *2,0 “1,8 


1 


can now be written 


a3 


<N,>{p; }da dp 


- p N, 0. 0 
<v;> 5 | SAjAQAz> [9 Nz, Of} <N;>{p;} da dB (6) 
=i = oi Se 


B. PROGRAM STRUCTURE 
The consistent pressure load program is structured as 
follows: 
Subroutine TRANS 
- reads element number, coordinates, and connectivity 
- reads the face to which pressure applied (NFACE) 
- reads pressure at each node of NFACE 
- reads number of Gauss points to be used 
in surface integration (2-6) 
- calls CUBAT 
Subroutine CUBAT 
- establishes vector of weights and sampling 
points to be used in Gaussian integration 
- calls SHAPE 
Subroutine SHAPE 
- evaluates shape functions at Gauss points 
- evaluates derivative of shape function at 
Gauss point 
- forms Jacobian matrix 
- evaluates determinant of Jacobian matrix to check 
for zero or negative Jacobian determinant 


- calls FACEP 
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Subroutine FACEP 
- establishes matrices in integrand 
- performs matrix multiplication to establish 
integrand 
- performs surface integration 
A listing of the consistent pressure load program is in 


Appendix C. 


C. PROGRAM VERIFICATION 

Several test cases were run to insure the accuracy of 
the pressure load preprocessor. A uniform distributed sur- 
face load was rotated through all faces of a one inch cube, 
and the 2-6 Gauss point integration scheme was checked for 
each face. The results agreed exactly, both in magnitude 
and sign, with the consistent load allocation found in 
Ref. [14]. 

The program was successfully checked to insure that dis- 
tributed surface force could be applied to more than one 
face of the element simultaneously, and that pressures could 
be applied simultaneously to faces that share a common edge. 

The ceramic flexure specimen mesh (Case 4) was also 
used to verify an actual problem. Pressure loads were 
applied to various faces of the specimen, including adjoin- 
ing faces. In each case, the program produced correct re- 
sults. 


As a final test, a distributed load was also placed 


across the span between the supports. A sketch of the 


problem is shown in Fig. 34. The results were again com- 
pared with simple beam theory. The results were excellent 


as shown in Fig. 35. 


D. USER'S MANUAL 

Figure 36 shows a hexahedral element in the natural co- 
ordinate system. The node numbering convention is also 
shown in Fig. 36. Table V shows the face numbering conven- 
tion used in the pressure load preprocessor. Each face has 
four corner nodes and four midside nodes. The pressure dis- 
tribution acting on an element face is defined by specifying 
the nodal pressure intensities at the corner and mid-side 
nodes. For example, if Face 4 were loaded Pus Pes Poy P35 


and Pay must be specified in that order. 


Poor Pise Pig» 
Appendix B lists the input data and results for the sample 
bar problem, shown in Fig. 36. 


Ls Control Card <10n5) 


Columns Variable Meaning 
1-5 NUMNP Total number of node 


points in problem. 
6-10 NUMEL Total number of 
elements in problem. 
11-15 NGP Number of Gauss points 
to use in surface 
integration. (2<NGP<6) 
16-20 NCUR Number of ADINA load 


curve that output loads 


refer to. (usually,NCUR=1) 
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21-25 NPLOAD Number of element 
faces with distributed 
loads. 

26-30 IPRINT 1 l-program will just 
read and print ADINA 
input deck 

31-35 IPRINT 2 l-program will not 
print element-by- 
element solution. | 

36-40 IPRINT 3 l-program will not | 
punch output deck of 
consistent loads. 

: 41-45 IPRINT 4 l-program will omit 

printing input data. 
2. Nodal Pressures Card (215, 8F5.0) 

1-5 NPREL (J) element number 

6-10 NFACE (J) face to which pressure 
applied. (See Table V) 

11-50 PRESS (I,J) nodal pressure 
intensities input here: 
corner nodes first, 
then midside. (See 
Table V) Positive pres- 
sure is directed in 


Same sense as outward 


i unit normal. 


E. USE OF PRESSURE LOAD PREPROCESSOR ON CP/CMS 


The pressure load preprocessor is ideally suited for use 
on CP/CMS. A compiled version of the program (filetype=text) 
should be stored on a private disk or workspace. A data 
file, such as FTO4F002, should also be defined. The input 
data deck (Appendix B) can then be read into this file [15]. 

It is suggested that an executive type routine also be 
used: 

& TYPEOUT ERROR 

FILEDEF 05 DSK FILE FTO4FO02 (PERM) 
FILEDEF 06 PRT (PERM) 

FILEDEF 07 PUN (PERM) 

VSET RDYMSG OFF 

$ PRESSURE 

This defines your input/output files and executes the 
program. Output is directed to the line printer, and card 
punch. 

The pressure program utilizes approximately 300K BYTES 
of core. The extra core should be requested at LOGIN by: 

LannnnPtt.~ 300K 
where: 
mnnn = user number 


tt = terminal number 
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VII. USING PSAP1 AS A DISPLACEMENT POSTPROCESSOR 


PSAP1 is capable of postprocessing displacements. The 
results can be presented in two ways: 

1. The deformed structure can be plotted. In this case, 
scaled displacements are added to the original geometry to 
produce the deformed coordinates. The deformed structure is 
then plotted using this new geometry. 

2. The displacements can be represented as scaled vec- 
tors. The scaled displacement vectors are plotted at origi- 
nal geometry nodes. 

The purpose of this chapter is to explain the modifica- 
tion of the ADINA code to store the displacements, and how 


to use the PSAP1 postprocessing capabilities. 


A. MODIFICATION OF ADINA 


Reference [16] states that a version of SAP IV had been 
modified to punch a deck of displacement cards from the out- 
put data. When dealing with large quantities of data, how- 
ever, storage of this data on direct access storage devices 
is preferable. Therefore, the ADINA code was modified to 
store the output displacements on device 58, in addition to 
the usual printed output. 

The flag for saving the displacements is the integer 
number 1 in column 80 of the ADINA master control card [4]. 


The flag variable name is ITP 58, and is passed in common 


/EAST/. 
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The job control language (JCL) for ADINA must also be 
modified to define device 58: 


//FTS8FO01 DD DISP=OLD,UNIT=3330,VOL=SER=DISKO, 
he DSN=Snnnn.pppp 


where: 
nnnn is the user number 
pppp is the name assigned by the user to identify 
the data set. 
Since an unformatted write statement is used to store 
the data on device 58, a variable blocksize must be used 
when allocating disk space. Reference [17] explains this 


procedure. 


B. MODIFICATIONS OF PSAP1 

There are presently two subroutines in PSAP1 that can 
be used to input displacement data to PSAP1. 

1. Subroutine DATA 9 was added to PSAP1 by Losh [16]. 
This subroutine reads displacement data from punched cards. 
Reference [(16] discusses its use, and should be consulted. 

2. Subroutine DATA 5 reads displacement data directly 
from device 58. A JCL card must be added to the PSAP1 JCL 
to define device 58: | 

//GO.FT58001 DD DISP=OLD,UNIT=3330,VOL=SER=DISKO, 

// DSN=Snnnn.pppp 
C. USER'S MANUAL FOR DISPLACEMENT POSTPROCESSING 

The PSAP1l User's Manual [11] applies with the following 
modifications: 


a) Namelist Option 
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NUDISP 
NVDISP 
NWDISP 
KDATA 


NVALUS 


b) 


KDISP 


IDMAG 


DMAGS 


Variable, Value 


Namelist Pict 


1 
1 


nnn 


n 


Description 


direction displacements input 


“ 


y direction displacements input 

z direction displacements input 
subroutine DATA 5 used to read 
displacement data from device 58. 
nnn is the integer number of dis- 
placements set to be plotted. 
This is equal to the number of 


node points. 


n=l: plot of deformed structure 

n=3; displacements represented 
by vectors at notes. 

n=l: direct magnification of 
displacement data by DMAGS. 

n=2: scaling of displacement 
data to a maximum value 
of DMAGS. 

ris a floating point number 

which indicates the magnifica- 

tion of the displacements. (0.1 


has been found appropriate.) 


Figures 37 and 38 are examples of displacement postpro- 


cessing. 


The structure is the bending specimen loaded in 


four-point flexure. 


Figure 37 illustrates the deformed structure option of 


PSAP1. Figure 38 shows the displacements of element 15 


plotted as scaled vectors at the node points. 
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CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 


1. The 12-node brick is a stiffer model of a structure 


than the fully quadratic 20-node brick. 


2. The accuracy of the results obtained using a 12-node 


brick depends on geometry, loading, and the orientation of 


the element. 


3. The flexure problem showed that two orientations of 


the 12-node brick are stiffer than the third. Resulting dis- 


placements and stresses differed significantly, and in all 


cases were lower than the 20-node brick and theory. 


4. The airfotl study showed that the maximum principal 


stresses predicted in the root section could differ by as 


much as 14.9% depending on the basic element used. The 12- 


node brick under-predicted the 20-node results. The use of 


the 12-node brick could therefore produce an unsafe design. 


B. RECOMMENDATIONS FOR FUTURE WORK 


As is true with many vrojects, this thesis raised many 


new questions. The following areas of future work are 


Suggested: 


1. The ADINA code should be reviewed with an eye toward 


Streamlining it. A basic code using only two- and three- 


dimensional elements, and the linear isotropic material 


A reduced version of ADINA 


model would be a useful tool. 


cevid then be installed on the time-sharing system, or 
other stand-alone system. 

2. The Field program, developed at NPS, should be ex- 
amined and expanded if necessary. This program would serve 
as an excellent thermal preprocessor. It should be modi- 
fied to produce a nodal temperature tape that could be 
used in ADINA, or other code capable of using temperature 

inputs. 

3. Aerodynamic pressure loads should be examined. A 
program that could calculate the nodal pressure intensi- 
ties for an element would provide the necessary input in- 
formation for the distributed load preprocessor. 

4. Centrifugal, pressure and thermal loadings should 


be superposed. 
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SPECIFIC FUEL CONSUMPTION, LB/HP-HR 


T76 -— STATE-OF-THE-ART UNCOOLED 
METAL TURBINE AT 1840°F 


CYCLE PRESSURE RATIO 


PROGRAM GOAL 
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IN POWER, 10 PERCENT TURBINE INLET 


TEMPERATURE, °F 
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T76 WITH CERAMIC 
COMPONENTS 
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Figure 1. Engine Design Parameters 
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a. 4 Node Plane Element 
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b. 6 Node Plane Element: Orientation l. 
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c. 6 Node Plane Element: Orientation 2. 
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d. 8 Node Plane Element 


Figure 3. Two-Dimensional Plane Element 
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SCHEMATIC SHOWING SELF-ALIGNING CAPABILITY 


Figure 4. Schematic of 4-Point 
Flexure Rig 


Pekan” Loken 


R 
A Re 


Figure 5a. Bend Specimen Free~-Body Diagram 
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Figure 5b. Bend Specimen Shear Diagram 
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Figure 5c. Bend Specimen Moment Diagram 


a 


SEEN ET 


z 


0.125" 


k-_- 155" 


Figure 6b. x-z View of Half-Span with Guided 
End Condition 
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Figure 6c. y-z View of Specimen with Guided 
Midplane Boundary Condition 
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PLANE OF 
SYMMETRY 


Figure 7. Cross-Sectional View of Deformed 
Specimen Showing x-z Plane of 
Symmetry 
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Figure 8a. 12-Node Brick Orientation for Case-1 
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Figure 8b. 12-Node Brick Orientation for Case~2 
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Figure 8c. 12-Node Brick Orientation for Case-3 
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Figure 8d. 20-Node Brick; Case-4 
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Figure 9a. Free Body Diagram for Flexure Case 
Study 
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Figure 9b. Free Body Diagram for Temperature Gradient 
Case Study 
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Figure 20. Trailing Edge Element Showing 
Re-entrant Angle 
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Figure 21. 12-Node Airfoil Section Showing 
Truncation of Leading Edge Element 
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TABLE I: MESH CHARACTERISTICS 


20-NODE BRICKS 
ElGQMGVNt CYP i ccc see we wesenctseedscceceencoenode brick 
Number Of nodeSerccccccccscccccccccccvestoa 
Number of eleMentsS..ccccescccrccccesceeeedd 
Degrees of fFreedom..cecccscccccccceceveellS2 
Mean half=bandwidth...cccescccccsccseeeeed3 
Maximum half=bandwidth....sccccccccsccceel29 


Number of stiffness matrix elements.....107064% 


12-NODE BRICKS 
ELGMGNE CyPGns cies veceessicieecdsvcsvewslL2=node brick 
Number of nodeSccccccccccccvecccesceseserl0 
Number of eleMentS..ccecesevescccesceseseda 
Degrees of Freedom. .cccceevececsccveee ed 48 
Mean half-bandwidth...cccccccsccccccvevetd 


Maximum half=-bandwidthececcccscccccceese0O 


Number of stiffness matrix elements.....29106 


TABLE II. COMPARISON OF MAXIMUM STRESSES IN AIRFOIL ROOT 
SECTION (Cf. Fig. 27) 
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TABLE III. COMPARISON OF MAXIMUM PRINCIPAL STRESSES IN 
AIRFOIL ROOT SECTION (cf. Fig. 27) 


al 
LOCATION |12-NODE AIRFOIL |20-NODE AIRFOIL|% DIFFERENCE 


9, (20) ~ 0, (12) 
% Difference = ——— 
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TABLE IV. COMPARISON OF SELECTED THERMAL STRESSES IN 
AIRFOIL ROOT SECTION (cf. Fig. 28) 


LOCATION 


% Difference = ote = _9(12) 
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APPENDIX A: PREPARATION OF TEMPERATURE TAPE 


If temperature input is used with the ADINA program, 
nodal point temperatures are stored in sequence on device 
number 56 [4]. At NPS, ITEL 3330 disks have been used, 
although any storage device may be designated. NSTE + 1 
records are read, where NSTE is the number of solution 
steps. (In linear static analysis, NSTE usually equals 1.) 

The read statement in the program is: 

READ (56) TIME, (TEMPV(I), I=1, NUMNP). 

On each record, the problem time of the nodal tempera- 
tures is read, followed by an implicit "do-loop" from 1 to 
the total number of node points (NUMNP) to read each nodal 
point temperature. 

A program to establish the temperature tape is listed 
at the end of this Appendix. It will be discussed in three 
parts: (A) job control language, (B) the RDADIN subroutine 


and (C) the binary write. 


A. TEMPERATURE TAPE JOB CONTROL LANGUAGE (JCL) 

The program listing in Appendix B starts with a two step 
procedure [(17]. The first is a purge operation that returns 
the previously reserved disk space to the operating system, 
the second step reserves disk space for the data set which 
is named FO099.TAPE. One cylinder of 3330 disk space (248K 


bytes) is reserved on one of the NPS computer center resident 


; 


PBK shite cae 


ans "cD aa wid lh 


tie ng tn alent nae 


disks. The data set has a variable block size, and an ex- 
piration date of 360 days after establishment. 
The JCL employed combines ease of continued use of the 


disk space, with minimum chances of data set expiration. 


B. THE RDADIN SUBROUTINE 

The RDADIN subroutine was originally programmed as part 
of the PSAP1 preprocessor. The subroutine will read an 
ADINA input deck from title card through element connectivity. 
It generates node numbers, nodal coordinates, element numbers 
and element connectivities for the entire mesh. Each of the 
above is stored in a linear array. The programmer is then 
able to utilize any or all of these arrays to input the 
nodal point temperature. 

For example, in the linear airfoil mesh there were iso- 
thermal planes. The matrix of z-coordinates (radial distance 
from axis of rotation) was used. Each z-coordinate was 
checked by a series of ‘if statements', and nodal tempera-~ 


tures were assigned. 


C. THE BINARY WRITE 

After the matrix of nodal temperatures was established, 
the next step was to write it sequentially onto device 56. 
Array RCD1 contained the problem time for the temperatures, 
and temperatures for each node; i.e., it was a linear array 
of dimension NUMNP + 1. This array could then be written 
onto device 56 by the statement: 


WRITE (56)RCD1. 


Device 56 must then be defined in the GO step by: 


//GO.FTS6F001 DD UNIT=3330,VOL=SER-DISKO1, 
//.DSN=F0099.TAPE,DISP=OLD 


Each time the WRITE statement is executed one record is 
established. The sequence must be repeated until NSTE + 1 
records are established, where NSTE is the number of solu- 
tion steps. 

In the linear airfoil problem, there was only one solu- 
tion time step. Therefore, only two records on device 56 
were needed. One record was at time=0, and the next record 
was one time step increment later. 

It is recommended that a hard copy of the RCD1 matrix 


be printed to ascertain the accuracy of the temperature 


tape. 
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APPENDIX B. NODAL TEMPERATURE TAPE PROGRAM 
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APPENDIX C. LISTING OF INPUT DECK AND RESULTS FOR 
SAMPLE PRESSURE PROBLEM 
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